Energy minimization for classification, pattern recognition, sensor fusion, data compression, network reconstruction and signal processing

ABSTRACT

A data analyzer/classifier comprises using a preprocessing step, and energy minimization step, and a postprocessing step to analyze/classify data.

REFERENCE TO RELATED APPLICATIONS

The present patent document is a a division of application Ser. No.11/097,783, filed Apr. 1, 2005, issued as U.S. Pat. No. 7,702,155 onApr. 20, 2010, which is a division of application Ser. No. 09/581,949,filed Jun. 19, 2000, issued as U.S. Pat. No. 6,993,186 on Jan. 31, 2006,which claims the benefit of the filing date under 35 U.S.C. §119(e) ofProvisional U.S. Patent Application Ser. No. 60/071,592, filed Dec. 29,1997. All of the foregoing applications and patents are herebyincorporated by reference.

COPYRIGHT NOTICE

A portion of the disclosure of this patent document contains materialwhich is subject to copyright protection. The copyright owner has noobjection to the facsimile reproduction by anyone of the patent documentor the patent disclosure, as it appears in the Patent and TrademarkOffice patent file or records, but otherwise reserves all copyrightrights whatsoever.

REFERENCE TO APPENDIX [CD ROM/SEQUENCE LISTING]

A computer program listing appendix is included containing computerprogram code listings on a CD-Rom pursuant to 37 C.F.R. 1.52(e) and ishereby incorporated herein by reference in its entirety. The totalnumber of compact discs is 1 (two duplicate copies are filed herewith).Each compact disc includes one (1) file entitled SOURCE CODE APPENDIX.Each compact disc includes 23,167 bytes. The creation date of thecompact disc is Mar. 24, 2005.

BACKGROUND

The present invention relates to recognition, analysis, andclassification of patterns in data from real world sources, events andprocesses. Patterns exist throughout the real world. Patterns also existin the data used to represent or convey or store information about realworld objects or events or processes. As information systems processmore real world data, there are mounting requirements to build moresophisticated, capable and reliable pattern recognition systems.

Existing pattern recognition systems include statistical, syntactic andneural systems. Each of these systems has certain strengths which lendsit to specific applications. Each of these systems has problems whichlimit its effectiveness.

Some real world patterns are purely statistical in nature. Statisticaland probabilistic pattern recognition works by expecting data to exhibitstatistical patterns. Pattern recognition by this method alone islimited. Statistical pattern recognizers cannot see beyond the expectedstatistical pattern. Only the expected statistical pattern can bedetected.

Syntactic pattern recognizers function by expecting data to exhibitstructure. While syntactic pattern recognizers are an improvement overstatistical pattern recognizers, perception is still narrow and thesystem cannot perceive beyond the expected structures. While some realworld patterns are structural in nature, the extraction of structure isunreliable.

Pattern recognition systems that rely upon neural pattern recognizersare an improvement over statistical and syntactic recognizers. Neuralrecognizers operate by storing training patterns as synaptic weights.Later stimulation retrieves these patterns and classifies the data.However, the fixed structure of neural pattern recognizers limits theirscope of recognition. While a neural system can learn on its own, it canonly find the patterns that its fixed structure allows it to see. Thedifficulties with this fixed structure are illustrated by the well-knownproblem that the number of hidden layers in a neural network stronglyaffects its ability to learn and generalize. Additionally, neuralpattern recognition results are often not reproducible. Neural nets arealso sensitive to training order, often require redundant data fortraining, can be slow learners and sometimes never learn. Mostimportantly, as with statistical and syntactic pattern recognitionsystems, neural pattern recognition systems are incapable of discoveringtruly new knowledge.

Accordingly, there is a need for an improved method and apparatus forpattern recognition, analysis, and classification which is notencumbered by preconceptions about data or models.

BRIEF SUMMARY

By way of illustration only, an analyzer/classifier process for datacomprises using energy minimization with one or more input matrices. Thedata to be analyzed/classified is processed by an energy minimizationtechnique such as individual differences multidimensional scaling(IDMDS) to produce at least a rate of change of stress/energy. Using therate of change of stress/energy and possibly other IDMDS output, thedata are analyzed or classified through patterns recognized within thedata. The foregoing discussion of one embodiment has been presented onlyby way of introduction. Nothing in this section should be taken as alimitation on the following claims, which define the scope of theinvention.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a diagram illustrating components of an analyzer according tothe first embodiment of the invention; and

FIG. 2 through FIG. 10 relate to examples illustrating use of anembodiment of the invention for data classification, patternrecognition, and signal processing.

DETAILED DESCRIPTION OF THE PRESENTLY PREFERRED EMBODIMENTS

The method and apparatus in accordance with the present inventionprovide an analysis tool with many applications. This tool can be usedfor data classification, pattern recognition, signal processing, sensorfusion, data compression, network reconstruction, and many otherpurposes. The invention relates to a general method for data analysisbased on energy minimization and least energy deformations. Theinvention uses energy minimization principles to analyze one to manydata sets. As used herein, energy is a convenient descriptor forconcepts which are handled similarly mathematically. Generally, thephysical concept of energy is not intended by use of this term but themore general mathematical concept. Within multiple data sets, individualdata sets are characterized by their deformation under least energymerging. This is a contextual characterization which allows theinvention to exhibit integrated unsupervised learning andgeneralization. A number of methods for producing energy minimizationand least energy merging and extraction of deformation information havebeen identified; these include, the finite element method (FEM),simulated annealing, and individual differences multidimensional scaling(IDMDS). The presently preferred embodiment of the invention utilizesindividual differences multidimensional scaling (IDMDS).

Multidimensional scaling (MDS) is a class of automated, numericaltechniques for converting proximity data into geometric data. IDMDS is ageneralization of MDS, which converts multiple sources of proximity datainto a common geometric configuration space, called the common space,and an associated vector space called the source space. Elements of thesource space encode deformations of the common space specific to eachsource of proximity data. MDS and IDMDS were developed for psychometricresearch, but are now standard tools in many statistical softwarepackages. MDS and IDMDS are often described as data visualizationtechniques. This description emphasizes only one aspect of thesealgorithms.

Broadly, the goal of MDS and IDMDS is to represent proximity data in alow dimensional metric space. This has been expressed mathematically byothers (see, for example, de Leeuw, J. and Heiser, W., “Theory ofmultidimensional scaling,” in P. R. Krishnaiah and L. N. Kanal, eds.,Handbook of Statistics, Vol. 2. North-Holland, N.Y., 1982) as follows.Let S be a nonempty finite set, p a real valued function on S×S,

p:S×S→R.

p is a measure of proximity between objects in S. Then the goal of MDSis to construct a mapping {tilde over (ƒ)} from S into a metric space(X, d),

{tilde over (ƒ)}:S→X,

such that p(i, j)=p_(ij)≅d({tilde over (ƒ)}(i), {tilde over (ƒ)}(i)),that is, such that the proximity of object i to object j in S isapproximated by the distance in X between {tilde over (ƒ)}(i) and {tildeover (ƒ)}(j). X is usually assumed to be n dimensional Euclidean spaceR^(n), with n sufficiently small.

IDMDS generalizes MDS by allowing multiple sources. For k=1, . . . , mlet S_(k) be a finite set with proximity measure p_(k), then IDMDSconstructs maps

{tilde over (ƒ)}_(k): S_(k)→X

such that p_(k)(i, j)=p_(ijk)≅d({tilde over (f)}_(k)(i), {tilde over(ƒ)}_(k)(j)), for k=1, . . . , m.

Intuitively, IDMDS is a method for representing many points of view. Thedifferent proximities p_(k) can be viewed as giving the proximityperceptions of different judges. IDMDS accommodates these differentpoints of view by finding different maps {tilde over (ƒ)}_(k) for eachjudge. These individual maps, or their image configurations, aredeformations of a common configuration space whose interpoint distancesrepresent the common or merged point of view.

MDS and IDMDS can equivalently be described in terms of transformationfunctions. Let P=(p_(ij)) be the matrix defined by the proximity p onS×S. Then MDS defines a transformation function

ƒ:p_(ij)→d_(ij)(X),

where d_(ij)(X)=d({tilde over (ƒ)}(i), {tilde over (ƒ)}(j)), with {tildeover (ƒ)} the mapping from S→X induced by the transformation function ƒ.Here, by abuse of notation, X={tilde over (ƒ)}(S), also denotes theimage of S under {tilde over (ƒ)}. The transformation function ƒ shouldbe optimal in the sense that the distances ƒ(p_(ij)) give the bestapproximation to the proximities p_(ij). This optimization criterion isdescribed in more detail below. IDMDS is similarly re-expressed; thesingle transformation ƒ is replaced by m transformations ƒ_(k). Note,these ƒ_(k) need not be distinct. In the following, the image of S_(k)under ƒ_(k) will be written X_(k).

MDS and IDMDS can be further broken down into so-called metric andnonmetric versions. In metric MDS or IDMDS, the transformations ƒ(ƒ_(k))are parametric functions of the proximities p_(ij)(p_(ijk)). NonmetricMDS or IDMDS generalizes the metric approach by allowing arbitraryadmissible transformations ƒ (ƒ_(k)), where admissible means theassociation between proximities and transformed proximities (also calleddisparities in this context) is weakly monotone:

p _(ij) <p _(kl) implies ƒ(p _(ij))≦ƒ(p _(kl)).

Beyond the metric-nonmetric distinction, algorithms for MDS and IDMDSare distinguished by their optimization criteria and numericaloptimization routines. One particularly elegant and publicly availableIDMDS algorithm is PROXSCAL See Commandeur, J. and Heiser, W.,“Mathematical derivations in the proximity scaling (PROXSCAL) ofsymmetric data matrices,” Tech. report no. RR-93-03, Department of DataTheory, Leiden University, Leiden, The Netherlands. PROXSCAL is a leastsquares, constrained majorization algorithm for IDMDS. We now summarizethis algorithm, following closely the above reference.

PROXSCAL is a least squares approach to IDMDS which minimizes theobjective function

${\sigma \left( {f_{1},\ldots \mspace{14mu},f_{m},X_{1},\ldots \mspace{14mu},X_{m}} \right)} = {\sum\limits_{k = 1}^{m}\; {\sum\limits_{i < j}^{n}\; {{w_{ijk}\left( {{f_{k}\left( p_{ijk} \right)} - {d_{ij}\left( X_{k} \right)}} \right)}^{2}.}}}$

σ is called the stress and measures the goodness-of-fit of theconfiguration distances d_(ij)(X_(k)) to the transformed proximitiesƒ_(k)(p_(ijk)). This is the most general form for the objectivefunction. MDS can be interpreted as an energy minimization process andstress can be interpreted as an energy functional. The w_(ijk) areproximity weights. For simplicity, it is assumed in what follows thatw_(ijk)=1 for all i, j, k.

The PROXSCAL majorization algorithm for MDS with transformations issummarized as follows.

1. Choose a (possibly random) initial configuration X⁰.

2. Find optimal transformations ƒ(p_(ij)) for fixed distancesd_(ij)(X⁰).

3. Compute the initial stress

${\sigma \left( {f,X^{0}} \right)} = {\sum\limits_{i < j}^{n}\; {\left( {{f\left( p_{ij} \right)} - {d_{ij}\left( X^{0} \right)}} \right)^{2}.}}$

4. Compute the Guttman transform X of X⁰ with the transformedproximities ƒ(p_(ij)). This is the majorization step.

5. Replace X⁰ with X and find optimal transformations ƒ(p_(ij)) forfixed distances d_(ij)(X⁰).

6. Compute σ(ƒ, X⁰).

7.Go to step 4 if the difference between the current and previous stressis not less than ε, some previously defined number. Stop otherwise.

For multiple sources of proximity data, restrictions are imposed on theconfigurations X_(k) associated to each source of proximity data in theform of the constraint equation X_(k)=ZW_(k).

This equation defines a common configuration space Z and diagonal weightmatrices W_(k). Z represents a merged or common version of the inputsources, while the W_(k) define the deformation of the common spacerequired to produce the individual configurations X_(k). The vectorsdefined by diag(W_(k)), the diagonal entries of the weight matricesW_(k), form the source space W associated to the common space Z.

The PROXSCAL constrained majorization algorithm for IDMDS withtransformations is summarized as follows. To simplify the discussion,so-called unconditional IDMDS is described. This means the mtransformation functions are the same: ƒ₁=f₂= . . . =ƒ_(m).

1. Choose constrained initial configurations X_(k) ⁰.

2. Find optimal transformations ƒ(p_(ijk)) for fixed distancesd_(ij)(X_(k) ⁰).

3. Compute the initial stress

${\sigma \left( {f,X_{1}^{0},\ldots \mspace{14mu},X_{m}^{0}} \right)} = {\sum\limits_{k = 1}^{m}\; {\sum\limits_{i < j}^{n}{\left( {{f\left( p_{ijk} \right)} - {d_{ij}\left( X_{k}^{0} \right)}} \right)^{2}.}}}$

4. Compute unconstrained updates X _(k) of X⁰ _(k) using the Guttmantransform with transformed proximities ƒ(p_(ijk)). This is theunconstrained majorization step.

5. Solve the metric projection problem by finding X_(k) ⁺ minimizing

${h\left( {X_{1},\ldots \mspace{14mu},X_{m}} \right)} = {\sum\limits_{k = 1}^{m}{t\; {r\left( {X_{k} - {\overset{\_}{X}}_{k}} \right)}^{l}\left( {X_{k} - {\overset{\_}{X}}_{k}} \right)}}$

subject to the constraints X_(k)=ZW_(k). This step constrains theupdated configurations from step 4.

6. Replace X_(k) ⁰ with X_(k) ⁺ and find optimal transformationsƒ(p_(ijk)) for fixed distances d_(ij)(X_(k) ⁰).

7. Compute σ(ƒ, X₁ ⁰, . . . , X_(m) ⁰).

8. Go to step 4 if the difference between the current and previousstress is not less than ε, some previously defined number. Stopotherwise.

Here, tr(A) and A′ denote, respectively, the trace and transpose ofmatrix A.

It should be pointed out that other IDMDS routines do not contain anexplicit constraint condition. For example, ALSCAL (see Takane, Y.,Young, F, and de Leeuw, J., “Nonmetric individual differencesmultidimensional scaling: an alternating least squares method withoptimal scaling features,” Psychometrika, Vol. 42, 1977) minimizes adifferent energy expression (stress) over transformations,configurations, and weighted Euclidean metrics. ALSCAL also producescommon and source spaces, but these spaces are computed throughalternating least squares without explicit use of constraints. Eitherform of IDMDS can be used in the present invention.

MDS and IDMDS have proven useful for many kinds of analyses. However, itis believed that prior utilizations of these techniques have notextended the use of these techniques to further possible uses for whichMDS and IDMDS have particular utility and provide exceptional results.Accordingly, one benefit of the present invention is to incorporate MDSor IDMDS as part of a platform in which aspects of these techniques areextended. A further benefit is to provide an analysis technique, part ofwhich uses IDMDS, that has utility as an analytic engine applicable toproblems in classification, pattern recognition, signal processing,sensor fusion, and data compression, as well as many other kinds of dataanalytic applications.

Referring now to FIG. 1, it illustrates an operational block diagram ofa data analysis/classifier tool 100. The least energy deformationanalyzer/classifier is a three-step process. Step 110 is a front end fordata transformation. Step 120 is a process step implementing energyminimization and deformation computations—in the presently preferredembodiment, this process step is implemented through the IDMDSalgorithm. Step 130 is a back end which interprets or decodes the outputof the process step 120. These three steps are illustrated in FIG. 1.

It is to be understood that the steps forming the tool 100 may beimplemented in a computer usable medium or in a computer system ascomputer executable software code. In such an embodiment, step 110 maybe configured as a code, step 120 may be configured as second code andstep 120 may be configured as third code, with each code comprising aplurality of machine readable steps or operations for performing thespecified operations. While step 110, step 120 and step 130 have beenshown as three separate elements, their functionality can be combinedand/or distributed. It is to be further understood that “medium” isintended to broadly include any suitable medium, including analog ordigital, hardware or software, now in use or developed in the future.

Step 110 of the tool 100 is the transformation of the data into matrixform. The only constraint on this transformation for the illustratedembodiment is that the resulting matrices be square. The type oftransformation used depends on the data to be processed and the goal ofthe analysis. In particular, it is not required that the matrices beproximity matrices in the traditional sense associated with IDMDS. Forexample, time series and other sequential data may be transformed intosource matrices through straight substitution into entries of symmetricmatrices of sufficient dimensionality (this transformation will bediscussed in more detail in an example below). Time series or othersignal processing data may also be Fourier or otherwise analyzed andthen transformed to matrix form.

Step 120 of the tool 100 implements energy minimization and extractionof deformation information through IDMDS. In the IDMDS embodiment of thetool 100, the stress function σ defines an energy functional overconfigurations and transformations. The configurations are furtherrestricted to those which satisfy the constraint equations X_(k)=ZW_(k).For each configuration X_(k), the weight vectors diag(W_(k)) are thecontextual signature, with respect to the common space, of the k-thinput source. Interpretation of σ as an energy functional isfundamental; it greatly expands the applicability of MDS as an energyminimization engine for data classification and analysis.

Step 130 consists of both visual and analytic methods for decoding andinterpreting the source space W from step 120. Unlike traditionalapplications of IDMDS, tool 100 often produces high dimensional output.Among other things, this makes visual interpretation and decoding of thesource space problematic. Possible analytic methods for understandingthe high dimensional spaces include, but are not limited to, linearprogramming techniques for hyperplane and decision surface estimation,cluster analysis techniques, and generalized gravitational modelcomputations. A source space dye-dropping or tracer technique has beendeveloped for both source space visualization and analyticpostprocessing. Step 130 may also consist in recording stress/energy, orthe rate of change of stress/energy, over multiple dimensions. The graphof energy (rate or change or stress/energy) against dimension can beused to determine network and dynamical system dimensionality. The graphof stress/energy against dimensionality is traditionally called a screenplot. The use and purpose of the scree plot is greatly extended in thepresent embodiment of the tool 100.

Let S={S_(k)} be a collection of data sets or sources S_(k) for k=1, . .. , m. Step 110 of the tool 100 converts each S_(k) ∈ S to matrix formM(S_(k)) where M(S_(k)) is a p dimensional real hollow symmetric matrix.Hollow means the diagonal entries of M(S_(k)) are zero. As indicatedabove, M(S_(k)) need not be symmetric or hollow, but for simplicity ofexposition these additional restrictions are adopted. Note also that thematrix dimensionality p is a function of the data S and the goal of theanalysis. Since M(S_(k)) is hollow symmetric, it can be interpreted andprocessed in IDMDS as a proximity (dissimilarity) matrix. Step 110 canbe represented by the map

M:S→H^(p)(R),

S_(k)→M(S_(k))

where )H^(p) (R is the set of p dimensional hollow real symmetricmatrices. The precise rule for computing M depends on the type of datain S, and the purpose of the analysis. For example, if S contains timeseries data, then M might entail the straightforward entry-wise encodingmentioned above. If S consists of optical character recognition data, orsome other kind of geometric data, then M(S_(k)) may be a standarddistance matrix whose ij-th entry is the Euclidean distance between “on”pixels i and j. M can also be combined with other transformations toform the composite, (M∘F)(S_(k)), where F, for example, is a fastFourier transform (FFT) on signal data S_(k). To make this moreconcrete, in the examples below M will be explicitly calculated in anumber of different ways. It should also be pointed out that for certaindata collections S it is possible to analyze the conjugate or transposeS′ of S. For instance, in data mining applications, it is useful totranspose records (clients) and fields (client attributes) thus allowinganalysis of attributes as well as clients. The mapping M is simplyapplied to the transposed data.

Step 120 of the presently preferred embodiment of the tool 100 is theapplication of IDMDS to the set of input matrices M(S)={M(S_(k))}. EachM(S_(k)))∈M(S) is an input source for IDMDS. As described above, theIDMDS output is a common space Z⊂R^(n) and a source space W. Thedimensionality n of these spaces depends on the input data S and thegoal of the analysis. For signal data, it is often useful to set n=p−1or even n=|S_(k)| where |S_(k)| denotes the cardinality of S_(k). Fordata compression, low dimensional output spaces are essential. In thecase of network reconstruction, system dimensionality is discovered bythe invention itself.

IDMDS can be thought of as a constrained energy minimization process. Asdiscussed above,the stress σ is an energy functional defined overtransformations and configurations in R^(n); the constraints are definedby the constraint equation X_(k)=ZW_(k). IDMDS attempts to find thelowest stress or energy configurations X_(k) which also satisfy theconstraint equation. (MDS is the special case when each W_(k)=I, theidentity matrix.) Configurations X_(k) most similar to the sourcematrices M(S_(k)) have the lowest energy. At the same time, each X_(k)is required to match the common space Z up to deformation defined by theweight matrices W_(k). The common space serves as a characteristic, orreference object. Differences between individual configurations areexpressed in terms of this characteristic object with these differencesencoded in the weight matrices W_(k). The deformation informationcontained in the weight matrices, or, equivalently, in the weightvectors defined by their diagonal entries, becomes the signature of theconfigurations X_(k) and hence the sources S_(k) (through M(S_(k)))).The source space may be thought of as a signature classification space.

The weight space signatures are contextual; they are defined withrespect to the reference object Z. The contextual nature of the sourcedeformation signature is fundamental. As the polygon classificationexample below will show, Z-contextuality of the signature allows thetool 100 to display integrated unsupervised machine learning andgeneralization. The analyzer/classifier learns seamlessly and invisibly.Z-contextuality also allows the tool 100 to operate without a prioridata models. The analyzer/classifier constructs its own model of thedata, the common space Z.

Step 130, the back end of the tool 100, decodes and interprets thesource or classification space output W from IDMDS. Since this outputcan be high dimensional, visualization techniques must be supplementedby analytic methods of interpretation. A dye-dropping or tracertechnique has been developed for both visual and analyticpostprocessing. This entails differential marking or coloring of sourcespace output. The specification of the dye-dropping is contingent uponthe data and overall analysis goals. For example, dye-dropping may betwo-color or binary allowing separating hyperplanes to be visually oranalytically determined. For an analytic approach to separatinghyperplanes using binary dye-dropping see Bosch, R. and Smith, J,“Separating hyperplanes and the authorship of the disputed federalistpapers,” American Mathematical Monthly, Vol. 105, 1998. Discretedye-dropping allows the definition of generalized gravitationalclustering measures of the form

${g_{p}\left( {A,x} \right)} = {\frac{\sum\limits_{y \neq x}\; {{\chi_{A}(x)}{\exp \left( {p \cdot {d\left( {x,y} \right)}} \right)}}}{\sum\limits_{y \neq x}{\exp \left( {p \cdot {d\left( {x,y} \right)}} \right)}}.}$

Here, A denotes a subset of W (indicated by dye-dropping), χ_(A)(x), isthe characteristic function on A, d(.,.) is a distance function, andp∈R. Such measures may be useful for estimating missing values in databases. Dye-dropping can be defined continuously, as well, producing akind of height function on W. This allows the definition of decisionsurfaces or volumetric discriminators. The source space W is alsoanalyzable using standard cluster analytic techniques. The preciseclustering metric depends on the specifications and conditions of theIDMDS analysis in question.

Finally, as mentioned earlier, the stress/energy and rate of change ofstress/energy can be used as postprocessing tools. Minima or kinks in aplot of energy, or the rate of change of energy, over dimension can beused to determine the dimensionality of complex networks and generaldynamical systems for which only partial output information isavailable. In fact, this technique allows dimensionality to be inferredoften from only a single data stream of time series of observed data.

A number of examples are presented below to illustrate the method andapparatus in accordance with the present invention. These examples areillustrative only and in no way limit the scope of the method orapparatus.

Example A Classification of Regular Polygons

The goal of this experiment was to classify a set of regular polygons.The collection S={S₁, . . . , S₁₆} with data sets S₁-S₄, equilateraltriangles; S₅-S₈, squares; S₉-S₁₂, pentagons; and S₁₃-S₁₆; hexagons.Within each subset of distinct polygons, the size of the figures isincreasing with the subscript. The perimeter of each polygon S_(k), wasdivided into 60 equal segments with the segment endpoints orderedclockwise from a fixed initial endpoint. A turtle application was thenapplied to each polygon to compute the Euclidean distance from eachsegment endpoint to every other segment endpoint (initial endpointincluded). Let x_(s) _(k) ^(i) denote the i-th endpoint of polygonS_(k), then the mapping M is defined by

M:S→H⁶⁰(R),

S_(k)→[d_(S) _(k) ₁ |d_(s) _(k) ₂ | . . . |d_(s) _(k) ₆₀ ]

where the columns

d _(s) _(k) _(i) =(d(x _(s) _(k) ^(i) ,x _(s) _(k) ¹),d(x _(s) _(k) ^(i),x _(s) _(k) ²), . . . ,d(x _(s) _(k) ^(i) ,x _(s) _(k) ⁶⁰))′.

The individual column vectors d_(s) _(k) _(i) have intrinsic interest.When plotted as functions of arc length they represent a geometricsignal which contains both frequency and spatial information.

The 16, 60×60 distance matrices were input into a publicly distributedversion of PROXSCAL. PROXSCAL was run with the following technicalspecifications: sources—16, objects—60, dimension—4, model—weighted,initial configuration—Torgerson, conditionality—unconditional,transformations—numerical, rate of convergence—0.0, number ofiterations—500, and minimum stress—0.0.

FIG. 2 and FIG. 3 show the four dimensional common and source spaceoutput. The common space configuration appears to be a multifacetedrepresentation of the original polygons. It forms a simple closed pathin four dimensions which, when viewed from different angles, or, what isessentially the same thing, when deformed by the weight matrices,produces a best, in the sense of minimal energy, representation of eachof the two dimensional polygonal figures. The most successful suchrepresentation appears to be that of the triangle projected onto theplane determined by dimensions 2 and 4.

In the source space, the different types of polygons are arranged, andhence, classified, along different radii. Magnitudes within each suchradial classification indicate polygon size or scale with the smallerpolygons located nearer the origin.

The contextual nature of the polygon classification is embodied in thecommon space configuration. Intuitively, this configuration looks like asingle, carefully bent wire loop. When viewed from different angles, asencoded by the source space vectors, this loop of wire looks variouslylike a triangle, a square, a pentagon, or a hexagon.

Example B Classification of Non-Regular Polygons

The polygons in Example A were regular. In this example, irregularpolygons S={S₁, . . . , S₆} are considered, where S₁-S₃ are trianglesand S₄-S₆ rectangles. The perimeter of each figure S_(k) was dividedinto 30 equal segments with the preprocessing transformation M computedas in Example A. This produced 6, 30×30 source matrices which were inputinto PROXSCAL with technical specifications the same as those aboveexcept for the number of sources, 6, and objects, 30.

FIG. 4 and FIG. 5 show the three dimensional common and source spaceoutputs. The common space configuration, again, has a “holographic” orfaceted quality; when illuminated from different angles, it representseach of the polygonal figures. As before, this change of viewpoint isencoded in the source space weight vectors. While the weight vectorsencoding triangles and rectangles are no longer radially arranged, theycan clearly be separated by a hyperplane and are thus accuratelyclassified by the analysis tool as presently embodied.

It is notable that two dimensional IDMDS outputs were not sufficient toclassify these polygons in the sense that source space separatinghyperplanes did not exist in two dimensions.

Example C Time Series Data

This example relates to signal processing and demonstrates the analysistool's invariance with respect to phase and frequency modification oftime series data. It also demonstrates an entry-wise approach tocomputing the preprocessing transformation M.

The set S={S₁, . . . , S₁₂} consisted of sine, square, and sawtoothwaveforms. Four versions of each waveform were included, each modifiedfor frequency and phase content. Indices 1-4 indicate sine, 5-8 square,and 9-12 sawtooth frequency and phase modified waveforms. All signalshad unit amplitude and were sampled at 32 equal intervals x, for 0≦x≦2π.

Each time series S_(k), was mapped into a symmetric matrix as follows.First, an “empty” nine dimensional, lower triangular matrixT_(k)=(t_(ij) ^(k))=T(S_(k)) was created. “Empty” meant that T_(k) hadno entries below the diagonal and zeros everywhere else. Nine dimensionswere chosen since nine is the smallest positive integer m satisfying theinequality m(m−1)/2≧32 and m(m−1)/2 is the number of entries below thediagonal in an m dimensional matrix. The empty entries in T_(k) werethen filled in, from upper left to lower right, column by column, byreading in the time series data from S_(k). Explicitly: s₁ ^(k)=t₂₁^(k), the first sample in S_(k) was written in the second row, firstcolumn of T_(k); s₂ ^(k)=t₃₁ ^(k), the second sample in S_(k) waswritten in the third row, first column of T_(k), and so on. Since therewere only 32 signal samples for 36 empty slots in T_(k) the fourremaining entries were designated missing by writing—2 in thesepositions (These entries are then ignored when calculating the stress).Finally, a hollow symmetric matrix was defined by setting

M(S _(k))=T _(k) +T _(k) ^(t).

This preprocessing produced 12, 9×9 source matrices which were input toPROXSCAL with the following technical specifications: sources—12,objects—9, dimension—8, model—weighted, initial configuration—Torgerson,conditionality—unconditional, transformations—ordinal, approach toties—secondary, rate of convergence—0.0, number of iterations—500, andminimum stress—0.0. Note that the data, while metric or numeric, wastransformed as if it were ordinal or nonmetric. The use of nonmetricIDMDS has been greatly extended in the present embodiment of the tool100.

FIG. 6 shows the eight dimensional source space output for the timeseries data. The projection in dimensions seven and eight, as detailedin FIG. 7, shows the input signals are separated by hyperplanes intosine, square, and sawtooth waveform classes independent of the frequencyor phase content of the signals.

Example D Sequences, Fibonacci, etc.

The data set S={S₁, . . . , S₉} in this example consisted of ninesequences with ten elements each; they are shown in Table 1, FIG. 8.Sequences 1-3 are constant, arithmetic, and Fibonacci sequencesrespectively. Sequences 4-6 are these same sequences with some error ornoise introduced. Sequences 7-9 are the same as 1-3, but the negativel′s indicate that these elements are missing or unknown.

The nine source matrices M(S_(k))=(m_(ij) ^(k)) were defined by

m _(ij) ^(k) =|s _(i) ^(k) −s _(j) ^(k)|,

the absolute value of the difference of the i-th and j-th elements insequence S_(k). The resulting 10×10 source matrices where input toPROXSCAL configured as follows: sources—9, objects—10, dimension—8,model—weighted, initial configuration—simplex,conditionality—unconditional, transformations—numerical, rate ofconvergence—0.0, number of iterations—500, and minimum stress—0.0.

FIG. 9 shows dimensions 5 and 6 of the eight dimensional source spaceoutput. The sequences are clustered, hence classified, according towhether they are constant, arithmetic, or Fibonacci based. Note that inthis projection, the constant sequence and the constant sequence withmissing element coincide, therefore only two versions of the constantsequence are visible. This result demonstrates that the tool 100 of thepresently preferred embodiment can function on noisy or errorcontaining, partially known, sequential data sets.

Example E Missing Value Estimation for Bridges

This example extends the previous result to demonstrate theapplicability of the analysis tool to missing value estimation on noisy,real-world data. The data set consisted of nine categories of bridgedata from the National Bridge Inventory (NBI) of the Federal HighwayAdministration. One of these categories, bridge material (steel orconcrete), was removed from the database. The goal was to repopulatethis missing category using the technique of the presently preferredembodiment to estimate the missing values.

One hundred bridges were arbitrarily chosen from the NBI. Each bridgedefined an eight dimensional vector of data with components the NBIcategories. These vectors were preprocessed as in Example D, creatingone hundred 8×8 source matrices. The matrices were submitted to PROXSCALwith specifications: sources—100, objects—8, dimension—7,model—weighted, initial configuration—simplex,conditionality—unconditional, transformations—numerical, rate ofconvergence—0.0, number of iterations—500, and minimum stress—0.00001.

The seven dimensional source space output was partially labeled bybridge material—an application of dye-dropping—and analyzed using thefollowing function

${g_{p}\left( {A_{i},x} \right)} = {\frac{\sum\limits_{y \neq x}\; {{\chi_{A_{i}}(x)} \cdot {d\left( {x,y} \right)}^{- p}}}{\sum\limits_{y \neq x}{d\left( {x,y} \right)}^{- p}}.}$

where p is an empirically determined negative number, d(x,y) isEuclidean distance on the source space, and χ₄ is the characteristicfunction on material set A_(i), i=1,2, where A₁ is steel, A₂ concrete.(For the bridge data, no two bridges had the same source spacecoordinates, hence g_(p) was well-defined.) A bridge was determined tobe steel (concrete) if g_(p)(A₁,x)>g_(p)(A₂,x)(g_(p)(A₁,x)<g_(p)(A₂,x)). The result was indeterminate in case ofequality.

The tool 100 illustrated in FIG. 1 estimated bridge constructionmaterial with 90 percent accuracy.

Example F Network Dimensionality for a 4-Node Network

This example demonstrates the use of stress/energy minima to determinenetwork dimensionality from partial network output data. Dimensionality,in this example, means the number of nodes in a network.

A four-node network was constructed as follows: generator nodes 1 to 3were defined by the sine functions, sin(2x), sin(2x+π/2), andsin(2x+4π/3); node 4 was the sum of nodes 1 through 3. The output ofnode 4 was sampled at 32 equal intervals between 0 and 2π.

The data from node 4 was preprocessed in the manner of Example D: theU-th entry of the source matrix for node 4 was defined to be theabsolute value of the difference between the i-th and j-th samples ofthe node 4 time series. A second, reference, source matrix was definedusing the same preprocessing technique, now applied to thirty two equalinterval samples of the function sin(x) for 0≦x≦2π. The resulting 2,32×32 source matrices were input to PROXSCAL with technicalspecification: sources—2, objects—32, dimension—1 to 6, model—weighted,initial configuration—simplex, conditionality—conditional,transformations—numerical, rate of convergence—0.0, number ofiterations—500, and minimum stress—0.0. The dimension specification hada range of values, 1 to 6. The dimension resulting in the loweststress/energy is the dimensionality of the underlying network.

Table 2, FIG. 10, shows dimension and corresponding stress/energy valuesfrom the analysis by the tool 100 of the 4-node network. Thestress/energy minimum is achieved in dimension 4, hence the tool 100 hascorrectly determined network dimensionality. Similar experiments wererun with more sophisticated dynamical systems and networks. Each ofthese experiments resulted in the successful determination of systemdegrees of freedom or dimensionality. These experiments included thedetermination of the dimensionality of a linear feedback shift register.These devices generate pseudo-random bit streams and are designed toconceal their dimensionality.

From the foregoing, it can be seen that the illustrated embodiment ofthe present invention provides a method and apparatus for classifyinginput data. Input data are received and formed into one or morematrices. The matrices are processed using IDMDS to produce astress/energy value, a rate or change of stress/energy value, a sourcespace and a common space. An output or back end process uses analyticalor visual methods to interpret the source space and the common space.The technique in accordance with the present invention therefore avoidslimitations associated with statistical pattern recognition techniques,which are limited to detecting only the expected statistical pattern,and syntactical pattern recognition techniques, which cannot perceivebeyond the expected structures. Further, the tool in accordance with thepresent invention is not limited to the fixed structure of neuralpattern recognizers. The technique in accordance with the presentinvention locates patterns in data without interference frompreconceptions of models or users about the data. The patternrecognition method in accordance with the present invention uses energyminimization to allow data to self-organize, causing structure toemerge. Furthermore, the technique in accordance with the presentinvention determines the dimension of dynamical systems from partialdata streams measured on those systems through calculation ofstress/energy or rate of change of stress/energy across dimensions.

While a particular embodiment of the present invention has been shownand described, modifications may be made. For example, PROXSCAL may bereplaced by other IDMDS routines which are commercially available or areproprietary. It is therefore intended in the appended claims to coverall such changes and modifications which fall within the true spirit andscope of the invention.

It is therefore intended that the foregoing detailed description beregarded as illustrative rather than limiting, and that it be understoodthat it is the following claims, including all equivalents, that areintended to define the spirit and scope of this invention.

1. A method for determining dimensionality of a network, thedimensionality corresponding to a number of degrees of freedom in thenetwork, the method comprisings: using a computer system, sampling datafrom one or more nodes of the network; using a computer system, mappingthe data into one or more matrices; using a computer system, applyingindividual differences multidimensional scaling to the one or morematrices to produce a common space output; and using a computer system,processing the common space output to determine the dimensionality ofthe network.
 2. The method for determining dimensionality of a networkof claim 1 wherein mapping the data into one or more matrices comprisesa processing step selected from the group consisting of forming one ormore symmetric matrices, forming one or more hollow symmetric matrices,entry-wise substituting the received data to populate one or moresymmetric matrices, forming one or more proximity matrices, forming oneor more distance matrices, and forming one or more Euclidean distancematrices.
 3. The method for determining dimensionality of a network ofclaim 1 wherein applying individual differences multidimensional scalingcomprises applying nonmetric individual differences multidimensionalscaling.
 4. A method for reconstructing a network, the methodcomprising: using a computer system, sampling data from one or morenodes of the network; using a computer system, mapping the data into oneor more matrices; using a computer system, applying individualdifferences multidimensional scaling to the one or more matrices toproduce a source space output; from the source space output, using acomputer system, determining the dimensionality of the network; using acomputer system, using free nodes to recreate and reconstruct individualnodes through the use of matrices containing missing values; and using acomputer system, establishing node connectivity through the use oflowest-energy connections constrained by dimensionality.
 5. The methodfor reconstructing a network of claim 4 wherein mapping the data intoone or more matrices comprises a processing step selected from the groupconsisting of forming one or more symmetric matrices, forming one ormore hollow symmetric matrices, entry-wise substituting the receiveddata to populate one or more symmetric matrices, forming one or moreproximity matrices, forming one or more distance matrices, and formingone or more Euclidean distance matrices.
 6. The method forreconstructing a network of claim 4 wherein applying applying individualdifferences multidimensional scaling comprises applying nonmetricindividual differences multidimensional scaling.
 7. A method fordetermining dimensionality of a dynamical system from partial data, thedimensionality corresponding to a number of degrees of freedom in thedynamical system, the method comprising: using a computer system,sampling data from the dynamical system; using a computer system,mapping the data into one or more matrices; using a computer system,applying individual differences multidimensional scaling to the one ormore matrices to produce a common space output; using a computer system,processing the common space output to determine dimensionality of thedynamical system.
 8. The method determining dimensionality of adynamical system from partial data of claim 7 wherein mapping the datainto one or more matrices comprises a processing step selected from thegroup consisting of forming one or more symmetric matrices, forming oneor more hollow symmetric matrices, entry-wise substituting the receiveddata to populate one or more symmetric matrices, forming one or moreproximity matrices, forming one or more distance matrices, and formingone or more Euclidean distance matrices.
 9. The method for determiningdimensionality of a dynamical system from partial data of claim 7wherein applying individual differences multidimensional scalingcomprises applying nonmetric individual differences multidimensionalscaling.